MODULE INIT
  USE GLOBVAR
  USE NRUTIL, only : ARRAY1D, NR_BIG2, NR_SMALL2
  use nr, only : sort
  USE INTERPOL, only : locate, ARRAY1D_INIT
  use random, only: set_seed, Sample_Normal, Sample_Uniform
  use MATRIX, only : Matrix_Inverse
  use STATISTICS, only : Mean
  use COMFUNCS
  
  use MPI  !For Fortran 90, working with Intel version of openmpi.
  IMPLICIT NONE
  
CONTAINS
  
  SUBROUTINE INIT_PARAS(id) !initialize parameter values which are independent of N_T
    IMPLICIT NONE
    INTEGER(4), INTENT(IN) :: id
    integer(4) :: i, age, ierr, j, iNend
    real(8) :: Vrtemp(19)
    
    !============================================================
    ! initialize other variables (mostly policy related)
    !============================================================        
    
    !The domain of parameters to be estimated
    !(pa_epsilon cannot be larger than 1 otherwise it dominates the value)
    !keep the upper bound of pa_epsilon at 0.6d0.
    !gr_paradomain(1,1:2) = (/0.0d0, 0.6d0/)   !pa_epsilon
    if (Ipt == 1) then
        gr_paradomain(1,1:2) = (/0.0d0, 0.6d0/)   !pa_epsilon, with part time option
    elseif (gcI_robust == 21 .OR. gcI_robust == 22 .OR. gcI_robust == 23) then
        gr_paradomain(1,1:2) = (/0.0d0, 1.0d0/)   !pa_epsilon  
    else
        if (gcL_searchBC) then
            gr_paradomain(1,1:2) = (/0.0d0, 1.0d0/)   !pa_epsilon
        else
            gr_paradomain(1,1:2) = (/0.0d0, 0.6d0/)   !pa_epsilon
        endif
    endif
    
    if (gcI_robust==12 .OR. gcI_robust==13) then  !exo Hone or lbd Hone
       gr_paradomain(2,1:2) = (/0.01d0, 0.14d0/)  !sigma (depreciation rate)
       gr_paradomain(3,1:2) = (/0.0d0, 1.0d0/)    !alpha_I, irrelevant
       gr_paradomain(4,1:2) = (/0.0d0, 1.0d0/)    !alpha_H, H (power on H)
    elseif (gcI_robust==8 .OR. gcI_robust==9 .OR. gcI_robust==10 .OR. gcI_robust==11) then  !exo or lbd
       gr_paradomain(2,1:2) = (/0.01d0, 0.14d0/)  !sigma (depreciation rate)
       gr_paradomain(3,1:2) = (/-0.1d0, 0.2d0/)    !alpha, I (coef on age)
       gr_paradomain(4,1:2) = (/-0.2d0, 0.01d0/)    !beta, H (coef on age^2)
    else
       if (gcI_robust == 21 .OR. gcI_robust == 22 .OR. gcI_robust == 23) then   !NoDep for all (21) or at work (22) or off work (23)
           gr_paradomain(2,1:2) = (/0.01d0, 0.14d0/)  !sigma (depreciation rate)
       else
           gr_paradomain(2,1:2) = (/0.01d0, 0.14d0/)  !sigma (depreciation rate)
       endif
       
       if (gcI_robust==5) then  !initial wealth = 50,000
          gr_paradomain(3,1:2) = (/0.01d0, 1.0d0/)    !alpha, on I
          gr_paradomain(4,1:2) = (/0.05d0, 1.0d0/)    !beta, on H 
       elseif (gcI_robust == 21 .OR. gcI_robust == 22 .OR. gcI_robust == 23) then  !no depreciation regardless of working or not (21); or if at work (22) or off work (23)
          gr_paradomain(3,1:2) = (/0.01d0, 0.99d0/)    !alpha, on I
          gr_paradomain(4,1:2) = (/0.0001d0, 1.0d0/)    !beta, on H 
       else
          gr_paradomain(3,1:2) = (/0.001d0, 1.0d0/)    !alpha, on I
          gr_paradomain(4,1:2) = (/0.001d0, 1.0d0/)    !beta, on H
       endif   
    endif
    i = 5 !for next parameter
    
    gr_paradomain(i,1:2) = (/1.0d-5, 2.5d0/)   !5 pXI_sd
    if (gcL_searchBC) then
        gr_paradomain(i+1,1:2) = (/-8.0d0, 5.0d0/) !6 (a0,pi) het mean of a0
    else
        gr_paradomain(i+1,1:2) = (/-10.0d0, -3.0d0/) !6 (a0,pi) het mean of a0
    endif
    gr_paradomain(i+2,1:2) = (/0.5d0, 3.0d0/)    !7 (a0,pi) het mean of pi
    i = i + 3  !8 for next parameter
    if ((N_het_a0==1) .and. (N_het_pi==1)) then
       gr_paradomain(i,1:2) = (/-10.0d0,10.0d0/) !8 H0_mean
       gr_paradomain(i+1,1:2) = (/0.0d0,5.0d0/) !9 H0_sd
       i = i + 2  !10 for next parameter
    else
       gr_paradomain(i,1:2) = (/1.0d-4, 5.0d0/)     !8 (a0,pi) het sd of a0
       gr_paradomain(i+1,1:2) = (/1.0d-4, 2.0d0/)   !9 (a0,pi) het sd of pi
       gr_paradomain(i+2,1:2) = (/-1.0d0, 1.0d0/)   !10 (a0,pi) het correlation rho
       gr_paradomain(i+3,1:2) = (/-1.0d0, 5.0d0/)   !11 H0_mean
       gr_paradomain(i+4,1:2) = (/0.0d0, 0.5d0/)    !12 H0_sd. was (/0.0d0, 3.0d0/)?
       i = i + 5  !13
       if (N_het_H0>1) then
          gr_paradomain(i,1:2) = (/-2.0d0, 5.0d0/)    !13 pinitH_a0
          gr_paradomain(i+1,1:2) = (/-2.0d0, 5.0d0/)  !14 pinitH_pi
          i = i + 2  !15
       endif
    endif
    if (gcL_searchBC) then
        gr_paradomain(i,1:2) = (/0.0d0, 5.0d6/)  !15 pb1
    else
        gr_paradomain(i,1:2) = (/1.0d3, 5.0d9/)  !15 pb1
    endif    
    i = i + 1       !i=16 for het case now
    if (N_pbequest==2) then
        gr_paradomain(i,1:2) = (/1.0d0, 1.0d3/)  !16 pb2
        i = i + 1       !i=17 for het case now
    endif
    
    if (gcL_searchBC) then
        gr_paradomain(i,1:2) = (/0.0d0, 1.0d0/)  !16 the ratio of borrowing constraint (0,1)
    else
        gr_paradomain(i,1:2) = (/0.1d0, 6.0d0/)  !16 petac !consumption, CRRA, 1.0d0 - (-3.0d0)
    endif
    
    !consumption shifter: 
    if (N_phi_t>0) then
        gr_paradomain(i+1,1:2) = (/-1.0d0, 5.0d0/)  !17 pphi1 (t)
        do j = 2, N_phi_t
            gr_paradomain(i+j,1:2) = (/-1.0d0, 1.0d0/)  !pphi2 (t^2): 18 (and) pphi3 (t^3): 19
        enddo !do j
    endif
    gr_paradomain(i+N_phi_t+1,1:2) = (/0.0d0, 5.0d0/)  !pphi(N_phi_t+1)  !on married. POSITIVE. 20
    
    !leisure taste, gamma_it
    gr_paradomain(i+N_phi_t+2,1:2) = (/-2.0d0, 5.0d0/)  !pa1  !taste of leisure for married, spouse not working 21
    gr_paradomain(i+N_phi_t+3,1:2) = (/-2.0d0, 5.0d0/)  !pa2  !taste of leisure, married, spouse working     22
    i = i + N_phi_t + 4   !next is 23
    
    if (IC_L_nonseparable >= 1) then
       gr_paradomain(i,1:2) = (/-1.23d0, 0.5d0/)  !pphi(N_phi_t+2): C_L nonseparable
       i = i + 1
    endif
    
    !extra utility for SSA between 62 and 70
    if (N_ssa_u>0) then
        gr_paradomain(i,1:2) = (/0.0d0, 0.1d0/)    !23, (/-0.1d0, 0.1d0/)  !pssa_62.   2.0E-4
        gr_paradomain(i+1,1:2) = (/0.0d0, 0.1d0/)  !24, (/-0.1d0, 0.1d0/)  !pssa_65.   2.0E-4
        gr_paradomain(i+2,1:2) = (/0.0d0, 0.1d0/)  !25, pssa_65_t.  5.0E-5 
        i = i + 3  !next is 26
    endif
    
    if (IHasHealthAge>=1) then
       ! the base for Good health 
       gr_paradomain(i,1:2) = (/-2.0d0, 0.0d0/)     !26, pa_health, taste of leisure for Excellent
      
       gr_paradomain(i+1,1:2) = (/0.0d0, 1.0d0/)   !27, pa_health, taste of leisure for bad health
       gr_paradomain(i+2,1:2) = (/0.0d0, 0.2d0/)    !28, pa_health_t 
       
       gr_paradomain(i+3,1:2) = (/1.0d0, 5.0d0/)   !29, pa_health, taste of leisure for disabled
       gr_paradomain(i+4,1:2) = (/0.0d0, 0.1d0/)    !30, pa_health_t
       i = i + 5  !index of next parameter
    endif
    
    if (Ipt == 1) then
       gr_paradomain(i,1:2) = (/-2.0d0, -0.5d0/)  !plpt_V(1), utility of working part time 26
       i = i + 1
       
       if (Ipt_coef==1 .OR. Ipt_coef==2 .OR. Ipt_coef==3 .OR. Ipt_coef==21 .OR. Ipt_coef==22) then
           gr_paradomain(i,1:2) = (/-0.1d0, 1.0d0/)  !plpt_V(2), t
           i = i + 1
           gr_paradomain(i,1:2) = (/-1.0d0, 1.0d0/)   !plpt_V(3), t^2 or alpha
           i = i + 1
       elseif (Ipt_coef==15) then !search center
           gr_paradomain(i,1:2) = (/-0.1d0, 1.0d0/)   !plpt_V(2), t
           i = i + 1
           gr_paradomain(i,1:2) = (/0.0d0, 100.0d0/)  !plpt_V(3), center
           i = i + 1
           gr_paradomain(i,1:2) = (/-1.0d0, 1.0d0/)   !plpt_V(4), t^2
           i = i + 1
       endif    
       
       if (Ipt_coef==21 .OR. Ipt_coef==22) then
           gr_paradomain(i,1:2) = (/-5.0d0, 5.0d0/)   !plpt_V(4), married spouse not working
           i = i + 1
           gr_paradomain(i,1:2) = (/-5.0d0, 5.0d0/)   !plpt_V(5), married spouse working
           i = i + 1 
       endif
    endif
    
    if (N_depreciation==2) then
        gr_paradomain(i,1:2) = (/0.01d0, 0.14d0/)  !psigma_offjob (depreciation rate off job)
        i = i + 1
    endif
    
    
    !Adult equivalence scale
    grAdultEquivScale = (/1.0d0, 1.34d0, 1.65d0, 1.97d0, 2.27d0/)
    
    !==========================================================================
    !  SIPP DATA - Moments. Only needed in the Master calculating GMM distance
    !==========================================================================
    if (id == master_mpi) then !MASTER
       MSM_Moments_Data%d1 = max(N_data_10, N_data_10_ss)
       MSM_Moments_Data%d2 = N_moments_type   !6+min(3,IHasHealthAge) +Ipt
       ALLOCATE( MSM_Moments_Data%a( MSM_Moments_Data%d1, MSM_Moments_Data%d2 ) )
       MSM_Moments_Data%a = -999.0d0  !-999 means missing.

       open(10,file=trim(cwd)//'../bpfcodes_com/data/raw_sipp_moments_age_v2.raw')
       j = 6
       
       do i = 1, 1000
          read(10,*,end=80) ierr, age, Vrtemp(1:j)   !1:lfpr_mean 2:lnw_mean 3:lnw_sd 4:lnw_fe_mean 5:E2U 6:U2E
          if (ierr .NE. Icol) cycle
          if ((age>=N_data_0) .AND. (age<=N_data_10)) then
             MSM_Moments_Data%a(age, 1:4) = Vrtemp(1:4)  !-999 means missing.
          elseif (age==1) then  !aggregate transition stored at moment=1, different from the data
             MSM_Moments_Data%a(1, 1) = Vrtemp(5)  !1.1. E to U
             MSM_Moments_Data%a(2, 1) = Vrtemp(6)  !2.1. U to E
          endif
       enddo
80     continue
       close(10)
       
       
       !5 consumption
       open(10,file=trim(cwd)//'../bpfcodes_com/data/raw_moments_consumption.raw')
       do i = 1, 1000
          read(10,*,end=81) age, Vrtemp(1:2)   !consumption for lowedu and highedu
          if ((age>=N_data_0_c) .AND. (age<=N_data_10)) then
             MSM_Moments_Data%a(age, 5) = Vrtemp(1+Icol)
          endif
       enddo
81     continue
       close(10)
       
       !6 ss
       open(10,file=trim(cwd)//'../bpfcodes_com/data/raw_sipp_moments_ss.raw')
       do i = 1, 1000
          read(10,*,end=82) ierr, age, Vrtemp(1:2)   !col age rSS frSS
          if (ierr .NE. Icol) cycle
          if ((age>=N_data_0_ss) .AND. (age<=N_data_10_ss)) then
             MSM_Moments_Data_ss(age, 1:2) = Vrtemp(1:2)   !col age rSS frSSp
             MSM_Moments_Data%a(age, 6) = Vrtemp(2)  !#6: frSSp (ssa)
          endif
       enddo
82     continue
       close(10)
       
       !3.1 wage depreciation after unemployment spell
       open(10,file=trim(cwd)//'../bpfcodes_com/data/raw_sipp_moments_unemp.raw')
       do i = 1, 10
          read(10,*,end=83) ierr, Vrtemp(1)   !col plnw41
          if (ierr .NE. Icol) cycle
          MSM_Moments_Data%a(3, 1) = Vrtemp(1)  !3.1 plnw41
       enddo
83     continue       
       
       !7 lfpr_diff_e2g
       !8 lfpr_diff_g2b
       !9 lfpr_diff_b2d
       if (IHasHealthAge>0) then
          j = 10 !for next moment 
          open(10,file=trim(cwd)//'../bpfcodes_com/data/raw_health_CPS_hs.raw')
          do i = 1, 1000
             read(10,*,end=85) age, Vrtemp(1:19)
             if ((age>=N_data_0_health) .AND. (age<=N_data_10)) then
                MSM_Moments_Data%a(age, 7) = Vrtemp(17)  !#7 LFPR difference Excellent-Good
                MSM_Moments_Data%a(age, 8) = Vrtemp(18)  !#8 LFPR difference Good-Bad
                MSM_Moments_Data%a(age, 9) = Vrtemp(19)  !#9 LFPR difference Bad-Disable
             endif
             !Health distribution and transition
             if ((age>=N_0) .AND. (age<=N_T)) then 
                gr_health_matrix(age,1,0) = Vrtemp(1)  !Distribution: excellent health
                gr_health_matrix(age,2,0) = Vrtemp(2)  !Distribution: good health
                gr_health_matrix(age,3,0) = Vrtemp(3)  !Distribution: bad health
                gr_health_matrix(age,4,0) = Vrtemp(4)  !Distribution: disab
                
                gr_health_matrix(age,1,1:4) = Vrtemp(5:8)  !transition: excellent to good/bad/disab
                gr_health_matrix(age,2,1:4) = Vrtemp(9:12)  !transition: good to good/bad/disab
                gr_health_matrix(age,3,1:4) = Vrtemp(13:16)  !transition: bad to good/bad/disab
                
                gr_health_matrix(age,4,1:3) = 0.0d0  !transition: disab is absorbing
                gr_health_matrix(age,4,4) = 1.0d0    !transition: disab is absorbing
             endif
          enddo
   85     continue
          close(10)
       else  !    IHasHealthAge==0
           j = 7 
       endif  !if (IHasHealthAge>0)   
       
       !part time working --- exclusive with 7 or 10 LFPR difference
       if (Ipt == 1) then
          open(10,file=trim(cwd)//'../bpfcodes_com/data/raw_sipp_moments_age_workpt.raw') !part time
          do i = 1, 1000
             read(10,*,end=86) ierr, age, Vrtemp(1)   !lfpr_pt
             if (ierr .NE. Icol) cycle
             if ((age>=N_data_0_pt) .AND. (age<=N_data_10_pt)) then
                MSM_Moments_Data%a(age, j) = Vrtemp(1)   ! 7 or 10
             endif
          enddo
   86     continue
          close(10)
       endif    !if (Ipt == 1)
       
       !==== MSM weight on GMM objective function =====
       MSM_wgt%d1 = N_moments_all
       MSM_wgt%d2 = MSM_wgt%d1
       allocate(MSM_wgt%a(MSM_wgt%d1, MSM_wgt%d2))
       !==== Variance-Covariance matrix of the data =====
       gmmse_S%d1 = MSM_wgt%d1
       gmmse_S%d2 = MSM_wgt%d1
       allocate(gmmse_S%a(gmmse_S%d1,gmmse_S%d2))
       gmmse_S%a = 0.0d0
       open(10,file=trim(cwd)//'../bpfcodes_com/data/raw_Cvarcov.raw')
       do i = 1, 3 * gmmse_S%d1 * gmmse_S%d2
          read(10,*, end=520) j, age, ierr, Vrtemp(1)
          if (j .NE. Icol) cycle
          if (Vrtemp(1) > -900.0d0) then
             gmmse_S%a(age,ierr) = Vrtemp(1)
             if (age .NE. ierr) gmmse_S%a(ierr,age) = Vrtemp(1)
          endif
       enddo
       520 continue
       close(10)
       
       !Moments: #5 consumption
       ierr = (N_data_10-N_data_0+1)*4
       do i = ierr+1, ierr+(N_data_10-N_data_0_c+1)
          gmmse_S%a(i, i) = 1.0d0 
       enddo  !do i
       
       !Moments: #6 SS Social Security
       ierr = ierr + (N_data_10-N_data_0_c+1)
       do i = ierr+1, ierr+(N_data_10_ss-N_data_0_ss+1)
          gmmse_S%a(i, i) = 1.0d0 
       enddo  !do i
    
       ierr = ierr + (N_data_10_ss-N_data_0_ss+1)
       !#7, #8, #9: lfpr_diff_e2g, lfpr_diff_g2b, lfpr_diff_b2d
       if (IHasHealthAge>0) then
           do i = ierr+1, ierr+(N_data_10-N_data_0_health+1)*3
               gmmse_S%a(i, i) = 1.0d0
           enddo  !do i
           ierr = ierr + (N_data_10-N_data_0_health+1)*3
       endif !if (IHasHealthAge>0)
       
       !Moments #7 or #10, part time
       if (Ipt == 1) then
           do i = ierr+1, ierr+(N_data_10_pt - N_data_0_pt +1)
               gmmse_S%a(i, i) = 1.0d0
           enddo  !do i
       endif
       
       !Moments: transition
       i = N_moments_all - 2 !E to U
       gmmse_S%a(i, i) = 0.01d0
       i = N_moments_all - 1 !U to E
       gmmse_S%a(i, i) = 0.05d0
       i = N_moments_all     !wage depreciation after unemployment
       gmmse_S%a(i, i) = 0.05d0
       
       ! OPTIMAL WEIGHT: only inverse of the diagonal.
       MSM_wgt%a = 0.0d0
       do i = 1, MSM_wgt%d1
          if (abs(gmmse_S%a(i,i)) > NR_EPS) then
             MSM_wgt%a(i,i) = 1.0d0 / gmmse_S%a(i,i)  
          else
             MSM_wgt%a(i,i) = 1.0d0 
          endif
       enddo
       
       !====================================================================================================
       !rescale the moments to have the same mean. This serves as an extra control for calculating distance
       !====================================================================================================
       MSM_wgt_scale%d1 = MSM_Moments_Data%d1  ! age
       MSM_wgt_scale%d2 = MSM_Moments_Data%d2  ! 1-6 or 9, or 1-7 or 10
       ALLOCATE(MSM_wgt_scale%a(MSM_wgt_scale%d1, MSM_wgt_scale%d2))
       MSM_wgt_scale%a = 1.0d0
       j = 0
       do i = 1, MSM_wgt_scale%d2  !1-6 or 9, 1-7 or 10
          iNend = N_data_10 
          if (i == 5) then  ! #5 consumption
             ierr = N_data_0_c
          elseif (i == 6) then ! #6 ss
             ierr = N_data_0_ss
             iNend = N_data_10_ss
          elseif ((IHasHealthAge>0) .AND. (i==7 .OR. i==8 .OR. i==9)) then ! #7, #8, #9 lfpr_diff
              ierr = N_data_0_health  !health
          elseif (Ipt==1 .AND. (i==7 .OR. i==10)) then
              ierr = N_data_0_pt
              iNend = N_data_10_pt
          else
             ierr = N_data_0  !22
          endif

          MSM_wgt_scale%a(ierr:iNend,i)=1.0d0/mean(MSM_Moments_Data%a(ierr:iNend,i)**2)

          do age = ierr, iNend
             select case (i)
                case (1)     !1 lfpr_mean
                   if (age<30+4*Icol) then
                      Vrtemp(1) = 200.0d0
                   elseif (age>50) then
                      Vrtemp(1) = 200.0d0
                   else
                      Vrtemp(1) = 200.0d0
                   endif
                case (2)     !2 lnw_mean
                   Vrtemp(1) = 2.5d0 * 200.0d0
                   if (age>50) Vrtemp(1) = Vrtemp(1) * 5.0d0
                case (3)     !3 lnw_sd
                   if (IHasHealthAge > 0) then
                       Vrtemp(1) = 1.0d0 
                   elseif (Icol == 1) then  !college
                       Vrtemp(1) = 1.0d0    
                   else
                       if (age<30) then
                           Vrtemp(1) = 1.0d0
                       else
                           Vrtemp(1) = 1.0d0
                       endif
                   endif  !if (IHasHealthAge > 0)
                case (4)     !4 lnw_fe_mean
                   Vrtemp(1) = 2.5d0 * 200.0d0
                   if (age>50) Vrtemp(1) = Vrtemp(1) * 5.0d0
                case (5)     !5 consumption.
                   if (gcI_robust==31 .OR. gcI_robust==32) then
                       Vrtemp(1) = 0.0d0
                   else
                       Vrtemp(1) = 10.0d0
                   endif    
                case (6)  !#6 ss
                    if (gcI_robust==31 .OR. gcI_robust==33) then
                        Vrtemp(1) = 0.0d0
                    else
                        Vrtemp(1) = 1.0d0
                    endif  
                case (7,8,9,10)     !#7, #8, #9 diff_lfpr, or diff_lfpr for Part Time
                   if (IHasHealthAge > 0 .AND. i>=7 .AND. i<=9) then
                      Vrtemp(1) = 1.0d0
                   elseif (Ipt==1 .and. (i==7 .OR. i==10)) then  !diff_lfpr for Part Time
                      Vrtemp(1) = 1.0d0 !Part Time
                   else
                      Vrtemp(1) = 0.0d0
                   endif
                case default
                   Vrtemp(1) = 0.0d0
            end select                
        
             j = j + 1
             !only re-scale the diagonal
             MSM_wgt%a(j,j) = MSM_wgt_scale%a(age,i) * Vrtemp(1)  !using MSM_wgt_scale only
          enddo !do age
       enddo  !i
       !aggregate transition rates
       j = N_moments_all - 2 !E to U
       MSM_wgt%a(j,j) = MSM_wgt%a(j,j) * 50.0d0
       j = N_moments_all - 1                        !U to E
       MSM_wgt%a(j,j) = MSM_wgt%a(j,j) * 50.0d0
       
       j = N_moments_all             ! wage depreciation after unemployment
       if (gcI_robust==8 .OR. gcI_robust==10 .OR. gcI_robust==12) then  !8: exo
           MSM_wgt%a(j,j) = 0.0d0
       elseif (gcI_robust==21) then !21:no H depreciation at all
           MSM_wgt%a(j,j) = 0.0d0
       elseif (IHasHealthAge > 0) then  !health
           MSM_wgt%a(j,j) = MSM_wgt%a(j,j) * 50.0d0    
       else
           MSM_wgt%a(j,j) = MSM_wgt%a(j,j) * 50.0d0
       endif
              
    endif  !if (id == master_mpi)
    
  END SUBROUTINE INIT_PARAS

 
 SUBROUTINE INIT_GRID(id,DimState_) !initialize grids - state, policy. Called only once in the beginning of the program
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: id, DimState_(:)
    integer(4) :: i1, i2, itemp, age, ierr, ipos(2), status(MPI_STATUS_SIZE)
    integer(2) :: iVageobs(60000,2)
    real(8) :: rminS, rmaxS, rtemp, Vrtemp(4)
 
    N_A = DimState_(1)    ! # of asset grids. large N_A is more important than N_H.
    N_H = DimState_(2)    ! # of H grids.
    N_AIME = DimState_(3) ! # of AIME grid.
    
    if (abs(grinitAsd)<1.0d-6) then
        VN_A(N_0) = 2      ! # of asset grids at N_0. Minimum 2
    else
        VN_A(N_0) = N_A
    endif
    
    VN_AIME(N_0) = 2   ! # of AIME grid at N_0. Minimum 2
    
    do i1 = N_0+1, N_T
        VN_A(i1) = min(N_A, max(N_A, VN_A(i1-1)+1))  ! # of asset grids. large N_A is more important than N_H.
        VN_AIME(i1) = min(N_AIME, max(N_AIME, VN_AIME(i1-1)+1)) ! # of AIME grid.
    enddo !do i1
    
    N_HEALTH = 1 + 3*min(1,IHasHealthAge)    !1 Excellent, 2 Good, 3 Bad, 4 Disable
    N_RECSS = 2

    !============================
    !  CPS: AIME share
    !============================
    cps_aimeshare%d1 = N_T
    allocate(cps_aimeshare%a( cps_aimeshare%d1 ))
    cps_aimeshare%a = 0.0d0
    
    !===========================================
    ! Initialize state spaces
    !===========================================
    ! Asset space.
    sgridA%d1 = N_T_grid
    sgridA%d2 = N_A
    allocate(sgridA%a(sgridA%d1, sgridA%d2))
    sgridA%a = 0.0d0
    
    sgridAnrlb%d1 = N_T_grid +1
    allocate(sgridAnrlb%a(sgridAnrlb%d1))
    sgridAnrlb%a = 0.0d0
    
    gr_Amax%d1 = sgridA%d1 !The maximum Asset level
    allocate(gr_Amax%a(gr_Amax%d1))
    gr_Amax%a = 0.0d0
    
    ! Human capital H space - H
    sgridH%d1 = sgridA%d1
    sgridH%d2 = N_H
    allocate(sgridH%a(sgridH%d1, sgridH%d2))
    sgridH%a = 0.0d0
    
    !AIME space, N_0:N_T
    sgridAIME%d1 = sgridA%d1
    sgridAIME%d2 = N_AIME
    allocate(sgridAIME%a(sgridAIME%d1, sgridAIME%d2))
    sgridAIME%a = 0.0d0
    
    !RECSS space
    sgridRECSS%d1 = N_RECSS
    ALLOCATE( sgridRECSS%a(sgridRECSS%d1) )
    sgridRECSS%a(1) = 0 !no ret, no ss 
    sgridRECSS%a(2) = 1 !ret, ss
    
       dvV_EV%d1 = N_A
       dvV_EV%d2 = N_AIME
       dvV_EV%d3 = N_HEALTH
       dvV_EV%d4 = N_RECSS
       dvV_EV%d5 = N_FS
       dvV_EV%d6 = N_H
       dvV_EV%d7 = N_T_grid
       ALLOCATE( dvV_EV%a( dvV_EV%d1, dvV_EV%d2, dvV_EV%d3, dvV_EV%d4, dvV_EV%d5, dvV_EV%d6, dvV_EV%d7) )
    
       do i1 = 0, Ipt
           !policy functions    
           dcI(i1)%d1 = N_A
           dcI(i1)%d2 = N_AIME
           dcI(i1)%d3 = N_HEALTH
           dcI(i1)%d4 = N_RECSS
           dcI(i1)%d5 = N_FS
           dcI(i1)%d6 = N_H
           dcI(i1)%d7 = dvV_EV%d7
           ALLOCATE( dcI(i1)%a( dcI(i1)%d1, dcI(i1)%d2, dcI(i1)%d3, dcI(i1)%d4, dcI(i1)%d5, dcI(i1)%d6, dcI(i1)%d7) )

           !threshold value
           depsilon(i1)%d1 = N_A
           depsilon(i1)%d2 = N_AIME
           depsilon(i1)%d3 = N_HEALTH
           depsilon(i1)%d4 = N_RECSS
           depsilon(i1)%d5 = N_FS
           depsilon(i1)%d6 = N_H
           depsilon(i1)%d7 = dvV_EV%d7
           ALLOCATE( depsilon(i1)%a( depsilon(i1)%d1, depsilon(i1)%d2, depsilon(i1)%d3,  &
                   & depsilon(i1)%d4, depsilon(i1)%d5, depsilon(i1)%d6, depsilon(i1)%d7) )
       enddo !do i1
       
       do i1 = 0, Ipt+1
          dvV(i1)%d1 = N_A
          dvV(i1)%d2 = N_AIME
          dvV(i1)%d3 = N_HEALTH
          dvV(i1)%d4 = N_RECSS
          dvV(i1)%d5 = N_FS
          dvV(i1)%d6 = N_H
          dvV(i1)%d7 = dvV_EV%d7
          ALLOCATE( dvV(i1)%a( dvV(i1)%d1, dvV(i1)%d2, dvV(i1)%d3, dvV(i1)%d4, dvV(i1)%d5, dvV(i1)%d6, dvV(i1)%d7) )
       
          dcc(i1)%d1 = N_A  !consumption policy function
          dcc(i1)%d2 = N_AIME
          dcc(i1)%d3 = N_HEALTH
          dcc(i1)%d4 = N_RECSS
          dcc(i1)%d5 = N_FS
          dcc(i1)%d6 = N_H
          dcc(i1)%d7 = dvV_EV%d7
          ALLOCATE( dcc(i1)%a( dcc(i1)%d1, dcc(i1)%d2, dcc(i1)%d3, dcc(i1)%d4, dcc(i1)%d5, dcc(i1)%d6, dcc(i1)%d7) )
       
          dciRECSS(i1)%d1 = N_A
          dciRECSS(i1)%d2 = N_AIME
          dciRECSS(i1)%d3 = N_HEALTH
          dciRECSS(i1)%d4 = N_RECSS
          dciRECSS(i1)%d5 = N_FS
          dciRECSS(i1)%d6 = N_H
          dciRECSS(i1)%d7 = dvV_EV%d7
          ALLOCATE(dciRECSS(i1)%a(dciRECSS(i1)%d1,dciRECSS(i1)%d2,dciRECSS(i1)%d3,dciRECSS(i1)%d4,dciRECSS(i1)%d5, &
                &  dciRECSS(i1)%d6,dciRECSS(i1)%d7))
          dciRECSS(i1)%a = 0       
       enddo
       
       !For passing values in each color, including master
       gI_HgridIncr = ceiling(dble(N_H) / dble(new_numprocs_mpi))  !used to distribute tasks among slaves and head slaves in the same group.
       gI_HgridFullNodes = N_H - new_numprocs_mpi * (gI_HgridIncr-1) !number of nodes calculating gI_HgridIncr points
    
       gI_VTasks_mpi%d1 = new_numprocs_mpi
       gI_VTasks_mpi%d2 = 4  !1-starting position on H, 2-end point, 3-count, 4-displacement
       allocate(gI_VTasks_mpi%a(gI_VTasks_mpi%d1,gI_VTasks_mpi%d2))
       gI_VTasks_mpi%a(1,1) = 1  !starting point
       gI_VTasks_mpi%a(1,3) = gI_HgridIncr !count
       gI_VTasks_mpi%a(1,2) = gI_VTasks_mpi%a(1,1)+gI_VTasks_mpi%a(1,3)-1 !end point
       gI_VTasks_mpi%a(1,4) = gI_displacement_mpi   !displacement
       do i1 = 2, new_numprocs_mpi
          gI_VTasks_mpi%a(i1,1) = gI_VTasks_mpi%a(i1-1,2) + 1
          if (i1<=gI_HgridFullNodes) then                
             gI_VTasks_mpi%a(i1,3) = gI_HgridIncr
          else
             gI_VTasks_mpi%a(i1,3) = gI_HgridIncr-1
          endif
          gI_VTasks_mpi%a(i1,2) = gI_VTasks_mpi%a(i1,1) + gI_VTasks_mpi%a(i1,3) - 1  !end point
          gI_VTasks_mpi%a(i1,4) = gI_VTasks_mpi%a(i1-1,4) + gI_VTasks_mpi%a(i1-1,3)  !displacement
       enddo !i
       gI_vsize_mpi = dcI(0)%d1*dcI(0)%d2*dcI(0)%d3*dcI(0)%d4*dcI(0)%d5
       
       if (myid_mpi==5) then
          open(10,file=trim(cwd)//'output/test_gI_VTasks_mpi.txt')
          do i1 = 1, gI_VTasks_mpi%d1
             write(10,'(I3,A,4I10)') i1, ", ", gI_VTasks_mpi%a(i1,:)
          enddo  !do i1
          close(10)
       endif   !if (myid_mpi==5)

    
    !=========================================================
    !  MSM: simulation
    !=========================================================
    MSM_Indiv%d1 = 22+2  !16
    MSM_Indiv%d2 = N_T_grid
    MSM_Indiv%d3 = N_Sim
    !1-A,2-H,3-AIME,4-c,5-I,6-L,7-Labor,8-Value,9-leisure error,10-XI error,
    !11-HEALTH shock,12-V0,13-V1, 14-Prob of leisure=0 or Phi(epsilon_tmin) in PT, 
    !15-FS error, 16-FS spouse inc
    !17-produced output, 18-ss tax, 19-medicare tax, 20-income tax, 21-received ss benefit, 22-ssdiinc
    !23-V(-1), 24-Phi(epsilon_tmax)
    ALLOCATE(MSM_Indiv%a(MSM_Indiv%d1, MSM_Indiv%d2, MSM_Indiv%d3))
    MSM_Indiv%a = 0.0d0
    MSM_Indiv_Irecss%d1 = N_T_grid
    MSM_Indiv_Irecss%d2 = N_Sim
    ALLOCATE(MSM_Indiv_Irecss%a(MSM_Indiv_Irecss%d1, MSM_Indiv_Irecss%d2))
    MSM_Indiv_Irecss%a = 0

    MSM_Indiv_Ihealth%d1 = N_T_grid
    MSM_Indiv_Ihealth%d2 = N_Sim
    ALLOCATE(MSM_Indiv_Ihealth%a(MSM_Indiv_Ihealth%d1, MSM_Indiv_Ihealth%d2))
    MSM_Indiv_Ihealth%a = 0
    
    MSM_Indiv_IFS%d1 = N_T_grid
    MSM_Indiv_IFS%d2 = N_Sim
    ALLOCATE(MSM_Indiv_IFS%a(MSM_Indiv_IFS%d1, MSM_Indiv_IFS%d2))
    MSM_Indiv_IFS%a = 0
    
    !Initial A (1) and H (2) level
    MSM_Indiv_initX%d1 = N_Sim
    MSM_Indiv_initX%d2 = 2  !1A, 2H
    ALLOCATE(MSM_Indiv_initX%a(MSM_Indiv_initX%d1, MSM_Indiv_initX%d2))
    MSM_Indiv_initX%a = 0.0d0
    !het type
    MSM_Indiv_Ihet%d1 = N_Sim
    MSM_Indiv_Ihet%d2 = 3 !!N_Sim; second dimension: 1. heterogeneity type, 2-observed first year, 3-observed last year
    ALLOCATE(MSM_Indiv_Ihet%a(MSM_Indiv_Ihet%d1, MSM_Indiv_Ihet%d2))
    
    ! map the het vectors, from (a0,pi) to group
    itemp = 0
    do i1 = 1, N_het_a0
        do i2 = 1, N_het_pi
            itemp = itemp + 1
            gI_het_pos(i1, i2) = itemp
            gI_het_group_p1(gI_het_pos(i1, i2)) = i1
            gI_het_group_p2(gI_het_pos(i1, i2)) = i2
        enddo !i2
    enddo !i1  

    !==== human capital space====
    sgridH%a = 0.0d0
    if (id == master_mpi) then !MASTER, importing H space data
       !Get the range, and then discretize for each period git.
       open(10,file=trim(cwd)//'../bpfcodes_com/data/raw_H_range.raw')  !18-100, rmin(18)=3.09
       do i1 = N_0, 200
          read(10,*,end=320) i2, rmaxS, rminS
          
          !was using 2.0d0*rmaxS, changed to 1.5
          if (i2>=N_0 .AND. i2<=N_T) sgridH%a(i2-N_0+1,1:2) = (/max(1.0d-2, rminS), 1.5d0*rmaxS/)
          
       enddo   
       320 continue
       close(10)
    endif !if (id == master_mpi)
    call MPI_BCAST(sgridH%a(1,1),sgridH%d1*2,MPI_DOUBLE_PRECISION,master_mpi,MPI_COMM_WORLD,ierr)
    sgridcom%d1 = N_H
    allocate(sgridcom%a(sgridcom%d1))
    call ARRAY1D_INIT(sgridcom,N_H,(/0.1d0,10.0d0/),1) !log linear
    rtemp = sgridcom%a(N_H)-sgridcom%a(1)
    do i1 = N_0-N_0+1, N_T-N_0+1
       
       if (i1==N_0-N_0+1) then
          rminS = sgridH%a(i1,1)
          rmaxS = sgridH%a(i1,2)
       else
          if (gcI_robust == 21) then   !21. depreciation is 0 for all
              rminS = sgridH%a(i1-1,1)
              rmaxS = max(3.0d0*sgridH%a(i1,2), sgridH%a(i1-1,2))
          else
              rminS = 0.86d0*sgridH%a(i1-1,1) !was using 0.85*
              if (gcI_robust == 8 .OR. gcI_robust == 10 .OR. gcI_robust == 12) then   !8. exo
                  rmaxS = max(3.0d0*sgridH%a(i1,2), sgridH%a(i1-1,2))
              elseif (gcI_robust == 21 .OR. gcI_robust == 22 .OR. gcI_robust == 23) then   !21 no dep at all, 22. no depreciation at work; 23 no depreciation off work;
                  rmaxS = max(2.0d0*sgridH%a(i1,2), sgridH%a(i1-1,2))    
              else
                  rmaxS = max(sgridH%a(i1,2), sgridH%a(i1-1,2))  !assume the H upperbound never decreases
              endif
          endif
       endif
       sgridH%a(i1,1) = rminS
       sgridH%a(i1,2) = rmaxS
    enddo
    do i1 = N_0-N_0+1, N_T-N_0+1
       rminS = sgridH%a(i1,1)
       rmaxS = sgridH%a(i1,2)
       if (Icol==1) rmaxS = rmaxS * 1.5d0 
       sgridH%a(i1,:) = rminS+(rmaxS-rminS)*(sgridcom%a-sgridcom%a(1))/rtemp
    enddo
    deallocate(sgridcom%a)     
    
    if (id == master_mpi) then !MASTER, importing other data
       !=== CPS: AIME share 
       i2 = 0
       open(10,file=trim(cwd)//'../bpfcodes_com/data/raw_aimeshare.raw')  !50-80
       do i1 = N_0, cps_aimeshare%d1
          read(10,*,end=100) age, rtemp
          if (age <= N_T) then
             cps_aimeshare%a(age) = rtemp
             if (i2<age) i2 = age
          endif
       enddo
       100 continue
       close(10)
       if (i2<=(N_T-2)) then
          cps_aimeshare%a(i2+1:N_T) = cps_aimeshare%a(i2) 
       elseif (i2==(N_T-1)) then
          cps_aimeshare%a(i2+1) = cps_aimeshare%a(i2)
       endif
       
       !=== The maximum Asset level, already normalized in the data
       open(10,file=trim(cwd)//'../bpfcodes_com/data/raw_A_max.raw') !18-100
       do i1 = N_0, 200
          read(10,*,end=310) i2, rmaxS
          if (i2 <= N_T) gr_Amax%a(i2-N_0+1) = max(6.0d0, rmaxS)
       enddo   
       310 continue
       close(10)
       
       !=== AIME space====
       !Get the range, and then discretize for each period git.
       open(10,file=trim(cwd)//'../bpfcodes_com/data/raw_AIME_range.raw')
       do i1 = N_0, 200
          read(10,*,end=330) i2, rmaxS, rminS
          rmaxS = max(rmaxS, 0.0d0)
          rminS = max(rminS, 0.0d0)
          if (i2>=N_0 .AND. i2<=N_T) then
             if (rmaxS < rminS) then
                rtemp = rmaxS
                rmaxS = rminS
                rminS = rtemp
             endif
             if (i2 == N_0) then
                rminS = 1.0d0/grTotHours  !set minimum to be 1.0d0/grTotHours, not 0.0d0
                rmaxS = 2.0d0 * rminS
             else
                rminS = sgridAIME%a(i2-1-N_0+1,1) !lower bound of AIME stays same across period
                if (gcI_robust == 21 .OR. gcI_robust == 22 .OR. gcI_robust == 23) then   !21 no dep at all, 22. no depreciation at work, 23 no depreciation off work
                    rtemp = 5.0d0
                else
                    rtemp = 3.0d0
                endif
                rmaxS = max(rtemp*rmaxS, sgridAIME%a(i2-1-N_0+1,2))
             endif
             sgridAIME%a(i2-N_0+1,1:2) = (/rminS, rmaxS/)
          endif
       enddo  !do i1 = N_0, 200
       330 continue
       close(10)
       
       do i1 = N_0, N_T
          if (VN_AIME(i1)==1) cycle  !sgridAIME%a(i1,1) is already rminS 
          
          sgridcom%d1 = VN_AIME(i1)
          allocate(sgridcom%a(sgridcom%d1))
          call ARRAY1D_INIT(sgridcom,VN_AIME(i1),(/0.1d0,10.0d0/),1) !0 linear; 1 log linear
          rtemp = sgridcom%a(VN_AIME(i1))-sgridcom%a(1) 
          rminS = sgridAIME%a(i1-N_0+1,1)
          rmaxS = sgridAIME%a(i1-N_0+1,2)
          if (Icol==1) rmaxS = rmaxS * (1.0d0 + 0.5d0*dble(i1-N_0)/dble(N_T-N_0))
          sgridAIME%a(i1-N_0+1,1:VN_AIME(i1)) = rminS+(rmaxS-rminS)*(sgridcom%a-sgridcom%a(1))/rtemp
          deallocate(sgridcom%a)
       enddo
    
       !Family status, transition and spouse income
       gr_FS_inc = 0.0d0
       open(10,file=trim(cwd)//'../bpfcodes_com/data/raw_sipp_spinc.raw')
       do i1 = 1, 1000
          read(10,*,end=340) ierr, age, Vrtemp(1:4)
          if (ierr .NE. Icol) cycle
          !1spsinc_bar 2spsinc_ln_bar 3sd_spsinc_bar 4sd_spsinc_ln_bar
          if ((age>=N_0) .AND. (age<=N_T)) then
             gr_FS_inc(age, 1) = Vrtemp(2) !mean of ln(inc)
             gr_FS_inc(age, 2) = Vrtemp(4) !sd of ln(inc)
          endif
       enddo
   340  continue
       close(10)
       !==== for XI Gauss-Hermite integration over spouse income
       call Gauss_Hermite(gr_FS_spinc_xw(N_0,:,1),gr_FS_spinc_xw(N_0,:,2)) !x,w
       !x = x * 1.41421356237d0 * sigma + mu
       gr_FS_spinc_xw(N_0,:,2) = gr_FS_spinc_xw(N_0,:,2) / 1.77245385091d0    ! /sqrt(pi)
       do i1 = N_0+1, N_T
          gr_FS_spinc_xw(i1,:,1) = gr_FS_spinc_xw(N_0,:,1)
          gr_FS_spinc_xw(i1,:,2) = gr_FS_spinc_xw(N_0,:,2)
       enddo
       do i1 = N_0, N_T
          gr_FS_spinc_xw(i1,:,3)=exp(gr_FS_spinc_xw(i1,:,1)*1.41421356237d0*gr_FS_inc(i1,2)+gr_FS_inc(i1,1)) !ln(spinc)
          do i2 = 1, N_FS_spinc_int
             gr_FS_spinc_xw(i1,i2,3) = max(0.0d0, gr_FS_spinc_xw(i1,i2,3))
          enddo !do i2
       enddo
       !1-x (standard normal), 2-weights, 3-actual spouse inc
       
       if (gcL_noFS) then
          do i1 = N_0, N_T
             gr_FS_spinc_xw(i1,:,3) = 0.0d0 
             do i2 = 1, N_FS 
                gr_FS_trans(i1, i2, 1) = 1.0d0
                gr_FS_trans(i1, i2, 2:) = 0.0d0
             enddo
          enddo
       else
          gr_FS_trans = 0.0d0
          open(10,file=trim(cwd)//'../bpfcodes_com/data/raw_sipp_Mt.raw')
          do i1 = 1, 1000
             read(10,*,end=350) ierr, age, i2, Vrtemp(1:3)
             if (ierr .NE. Icol) cycle
             if ((age>=N_0) .AND. (age<=N_T)) then
                gr_FS_trans(age, 1+i2, :) = Vrtemp(1:3) / sum(Vrtemp(1:3))
             endif
          enddo
   350   continue
          close(10) 
       endif  !if (gcL_noFS)
       
       
       !==== for XI Gauss-Hermite integration over XI
       call Gauss_Hermite(gr_XI_xw(:,1),gr_XI_xw(:,2)) !x,w
       !x = x * 1.41421356237d0 * sigma + mu
       gr_XI_xw(:,2) = gr_XI_xw(:,2) / 1.77245385091d0    ! /sqrt(pi)
       
       
       !==== het space: (a0,pi)
       call normal_discrete (N_het_a0, 0.0d0, 1.0d0, p1_Vhet_x(:,1), p1_Vhet_x(:,2))
       call normal_discrete (N_het_pi, 0.0d0, 1.0d0, p2_Vhet_x(:,1), p2_Vhet_x(:,2))
       
       !=== MSM: simulation ===
       ! This is in master.
       call Set_Seed(1234, 2345, 3456, 4567)
       
       do i1 = 1, N_Sim
          if (abs(grinitAsd)<1.0d-6) then
              MSM_Indiv_initX%a(i1, 1) = 0.0d0   !initial A level
          else
              MSM_Indiv_initX%a(i1, 1) = grinitAsd * Sample_Normal(0.0d0, 1.0d0)  !initial A level 
          endif
          
          MSM_Indiv_initX%a(i1, 2) = Sample_Normal(0.0d0, 1.0d0)  !initial H level
          !=============decide het type for each individual ===
          rminS = Sample_Uniform(0.0d0, 1.0d0)  !het: p1_x
          rtemp = 0.0d0
          do i2=1,N_het_a0-1
              rtemp = rtemp + p1_Vhet_x(i2,2)
              if (rminS<=rtemp) then
                  ipos(1) = i2
                  exit
              endif
          enddo !i2
          if ((N_het_a0 == 1) .OR. (rminS > rtemp)) ipos(1) = N_het_a0
          
          rminS = Sample_Uniform(0.0d0, 1.0d0)  !het: p2_x
          rtemp = 0.0d0
          do i2=1,N_het_pi-1
              rtemp = rtemp + p2_Vhet_x(i2,2)
              if (rminS<=rtemp) then
                  ipos(2) = i2
                  exit
              endif
          enddo !i2
          if ((N_het_pi == 1) .OR. (rminS > rtemp)) ipos(2) = N_het_pi
          
          !the MSM_Indiv_Ihet will be re-labeled later. Here it is just to count the number of each type
          MSM_Indiv_Ihet%a(i1,1) = gI_het_pos(ipos(1), ipos(2))
          
          !====================time varying errors
          do i2 = N_0, N_T-gcI_IncreaseLE
             MSM_Indiv%a(9,i2-N_0+1,i1) = Sample_Normal(0.0d0, 1.0d0)    !leisure error --- standard normal
             if (gcL_XI_discrete) then
                MSM_Indiv%a(10,i2-N_0+1,i1) = Sample_Uniform(0.0d0, 1.0d0)   !XI error  --- standard uniform, XI is discrete in simulation 
             else
                MSM_Indiv%a(10,i2-N_0+1,i1) = Sample_Normal(0.0d0, 1.0d0)   !XI error  --- standard normal
             endif
             MSM_Indiv%a(11,i2-N_0+1,i1) = Sample_Uniform(0.0d0, 1.0d0)   !HEALTH error  --- standard uniform
          enddo !do i2
       enddo !do i1
       
       
       !gI_Vcolor%a(,:) 1 number of procs, 2 number of simulations, 3 Incr, 4 # of Full nodes
       do i1 = 1, MSM_Indiv_Ihet%d1  !Here it is to count the number of each type
          gI_Vcolor%a(MSM_Indiv_Ihet%a(i1,1),2) = gI_Vcolor%a(MSM_Indiv_Ihet%a(i1,1),2) + 1
       enddo !do i1
       do i1 = 1, gI_Vcolor%d1
          gI_Vcolor%a(i1,3) = ceiling(dble(gI_Vcolor%a(i1,2)) / dble(gI_Vcolor%a(i1,1)))  !used to distribute tasks among slaves in the same group.
          gI_Vcolor%a(i1,4) = gI_Vcolor%a(i1,2) - gI_Vcolor%a(i1,1) * (gI_Vcolor%a(i1,3)-1) !number of nodes calculating gI_Vcolor%a(i1,3) points
       enddo
       !Re-label the MSM_Indiv_Ihet. Everything else is random, so it does not matter.
       itemp = 0
       do i1 = 1, gI_Vcolor%d1
          MSM_Indiv_Ihet%a(itemp+1:itemp+gI_Vcolor%a(i1,2),1) = i1
          itemp = itemp + gI_Vcolor%a(i1,2) 
       enddo
       
       if (gcI_IncreaseLE > 0) then  !this is to make sure shocks are same for the same years
          do i1 = 1, N_Sim
             do i2 = N_T-gcI_IncreaseLE+1, N_T
                MSM_Indiv%a(9,i2-N_0+1,i1) = Sample_Normal(0.0d0, 1.0d0)    !leisure error --- standard normal
                if (gcL_XI_discrete) then
                   MSM_Indiv%a(10,i2-N_0+1,i1) = Sample_Uniform(0.0d0, 1.0d0)   !XI error  --- standard uniform, XI is discrete in simulation 
                else
                   MSM_Indiv%a(10,i2-N_0+1,i1) = Sample_Normal(0.0d0, 1.0d0)   !XI error  --- standard normal
                endif
                MSM_Indiv%a(11,i2-N_0+1,i1) = Sample_Uniform(0.0d0, 1.0d0)   !HEALTH error  --- standard uniform
             enddo !do i2
          enddo !do i1
       endif
       
       !Family status
       do i1 = 1, N_Sim
          do i2 = N_0, N_T
             MSM_Indiv%a(15,i2-N_0+1,i1) = Sample_Uniform(0.0d0, 1.0d0) !FS transition error --- standard uniform
             if (gcL_spinc_discrete) then
                MSM_Indiv%a(16,i2-N_0+1,i1) = Sample_Uniform(0.0d0, 1.0d0)  !FS spouse income error --- standard uniform, discrete
             else
                MSM_Indiv%a(16,i2-N_0+1,i1) = Sample_Normal(0.0d0, 1.0d0)  !FS spouse income error --- standard normal, continuous
             endif
          enddo !do i2
       enddo !do i1
       
       
       itemp = 0
       open(10,file=trim(cwd)//'../bpfcodes_com/data/raw_sipp_ageobs.raw')
       do i1 = 1, 100000
          read(10,*,end=355) ierr, i2, ipos(1:2)
          if (ierr .NE. Icol) cycle
          if (ipos(2)-ipos(1)<2) cycle
          !1spsinc_bar 2spsinc_ln_bar 3sd_spsinc_bar 4sd_spsinc_ln_bar
          if ((ipos(1)>=N_data_0) .AND. (ipos(2)<=N_data_10)) then
              itemp = itemp + 1
              iVageobs(itemp,1:2) = ipos(1:2)
          endif
       enddo  !do i1
   355  continue
       close(10)
        
       !observed for a short period, age span >=2
       rtemp = dble(itemp)
       do i1 = 1, N_Sim
           i2 = ceiling(rtemp*Sample_Uniform(0.0d0,1.0d0))
           MSM_Indiv_Ihet%a(i1,2:3) = iVageobs(i2,1:2) 
           !N_Sim; second dimension: 1. heterogeneity type, 2-observed first year, 3-observed last year
       enddo !do i1
       
    endif  !if (id == master_mpi)

    ! These variables are for calculating the simulated moments, in all nodes
    !MSM: Simulated Moments from MSM
    MSM_Moments_Sim%d1 = N_T
    MSM_Moments_Sim%d2 = 6+min(3,IHasHealthAge) + Ipt  !6/7 or 9/10
    ALLOCATE(MSM_Moments_Sim%a(MSM_Moments_Sim%d1, MSM_Moments_Sim%d2))
        
    MSM_sim4m_i%d1 = N_T
    MSM_sim4m_i%d2 = MSM_Moments_Sim%d2
    allocate(MSM_sim4m_i%a(MSM_sim4m_i%d1,MSM_sim4m_i%d2))
        
    
    !!!!Master Broadcast to all head slaves and slaves
    !==== XI Gauss-Hermite integration
    call MPI_BCAST(gr_XI_xw,N_XI_int*3,MPI_DOUBLE_PRECISION,master_mpi,MPI_COMM_WORLD,ierr)
    !==== health transition
    call MPI_BCAST(gr_health_matrix(N_0,1,0),(N_T-N_0+1)*4*5,MPI_DOUBLE_PRECISION,master_mpi,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(p1_Vhet_x,N_het_a0*2,MPI_DOUBLE_PRECISION,master_mpi,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(p2_Vhet_x,N_het_pi*2,MPI_DOUBLE_PRECISION,master_mpi,MPI_COMM_WORLD,ierr)
    
    call MPI_BCAST(cps_aimeshare%a(1),cps_aimeshare%d1,MPI_DOUBLE_PRECISION,master_mpi,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(gr_Amax%a(1),gr_Amax%d1,MPI_DOUBLE_PRECISION,master_mpi,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(sgridAIME%a(1,1),sgridAIME%d1*sgridAIME%d2,MPI_DOUBLE_PRECISION,master_mpi,MPI_COMM_WORLD,ierr)
    
    call MPI_BCAST(gr_FS_inc,(N_T-N_0+1)*2,MPI_DOUBLE_PRECISION,master_mpi,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(gr_FS_trans,(N_T-N_0+1)*N_FS*N_FS,MPI_DOUBLE_PRECISION,master_mpi,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(gr_FS_spinc_xw,(N_T-N_0+1)*N_FS_spinc_int*3,MPI_DOUBLE_PRECISION,master_mpi,MPI_COMM_WORLD,ierr)
        
    call MPI_BCAST(MSM_Indiv_initX%a(1,1),MSM_Indiv_initX%d1*MSM_Indiv_initX%d2,MPI_DOUBLE_PRECISION,master_mpi,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(MSM_Indiv_Ihet%a(1,1),MSM_Indiv_Ihet%d1*MSM_Indiv_Ihet%d2,MPI_INTEGER,master_mpi,MPI_COMM_WORLD,ierr)
    call MPI_BCAST(MSM_Indiv%a(1,1,1),MSM_Indiv%d1*MSM_Indiv%d2*MSM_Indiv%d3,MPI_DOUBLE_PRECISION,master_mpi,MPI_COMM_WORLD,ierr)
    
    call MPI_BCAST(gI_Vcolor%a(1,1),gI_Vcolor%d1*gI_Vcolor%d2,MPI_INTEGER,master_mpi,MPI_COMM_WORLD,ierr)
    if (myid_mpi==5) then
       open(10,file=trim(cwd)//'output/test_Vcolor.txt')
       do i1 = 1, gI_Vcolor%d1
          write(10,'(A,I10,A,4I8)') "gI_Vcolor%a(", i1, ',1:4) = ', gI_Vcolor%a(i1,:)
       enddo
       close(10)
       
       open(10,file=trim(cwd)//'output/test_FS_spinc_xw.txt')
       do i1 = N_0, N_T
          write(10,'(I3,A,3F20.5)') i1, ", 1, ", gr_FS_spinc_xw(i1,:,1)
          write(10,'(I3,A,3F20.5)') i1, ", 2, ", gr_FS_spinc_xw(i1,:,2)
          write(10,'(I3,A,3F20.5)') i1, ", 3, ", gr_FS_spinc_xw(i1,:,3)
       enddo
       close(10)
       open(10,file=trim(cwd)//'output/test_FS_Mt.txt')
       do i1 = N_0, N_T
          write(10,'(I3,A,3F20.5)') i1, ", 1, ", gr_FS_trans(i1,1,:)
          write(10,'(I3,A,3F20.5)') i1, ", 2, ", gr_FS_trans(i1,2,:)
          write(10,'(I3,A,3F20.5)') i1, ", 3, ", gr_FS_trans(i1,3,:)
       enddo
       close(10)  
       open(10,file=trim(cwd)//'output/test_health_matrix.txt')
       do i1 = N_0, N_T
          write(10,'(I3,A,5F20.5)') i1, ", 1, ", gr_health_matrix(i1,1,:)
          write(10,'(I3,A,5F20.5)') i1, ", 2, ", gr_health_matrix(i1,2,:)
          write(10,'(I3,A,5F20.5)') i1, ", 3, ", gr_health_matrix(i1,3,:)
          write(10,'(I3,A,5F20.5)') i1, ", 4, ", gr_health_matrix(i1,4,:)
       enddo
       close(10)  
    endif
    
    ! This is for MPI_Allgatherv the individual simulations
    gI_Indiv_Tasks_mpi%d1 = numprocs_mpi  !including the master node
    gI_Indiv_Tasks_mpi%d2 = 4  !starting position on Indiv, end point, 3count, 4displayment
    allocate(gI_Indiv_Tasks_mpi%a(gI_Indiv_Tasks_mpi%d1,gI_Indiv_Tasks_mpi%d2)) 
    gI_Indiv_size_mpi = MSM_Indiv%d1*MSM_Indiv%d2
    
    !MSM_Indiv_Ihet
    
    !in all nodes
    !gI_Vcolor%a(,:) 1 number of procs, 2 number of simulations, 3 Incr, 4 # of Full nodes
    !gI_Indiv_Tasks_mpi%a(,:) 1 starting position on Indiv, 2 end point, 3count, 4displayment
    ipos(1) = 0  !overall id rank
    do i1 = 1, gI_Vcolor%d1  !number of colors 
       ipos(2) = 0  !id rank in the same color
       do i2 = 1, gI_Vcolor%a(i1,1) !number of nodes in this color
          ipos(1) = ipos(1) + 1
          ipos(2) = ipos(2) + 1
          if (ipos(1) == 1) then
             gI_Indiv_Tasks_mpi%a(ipos(1),1) = 1  !starting point
             gI_Indiv_Tasks_mpi%a(ipos(1),3) = gI_Vcolor%a(1,3) !count
             gI_Indiv_Tasks_mpi%a(ipos(1),2) = gI_Indiv_Tasks_mpi%a(ipos(1),1)+gI_Indiv_Tasks_mpi%a(ipos(1),3)-1 !end point
             gI_Indiv_Tasks_mpi%a(ipos(1),4) = gI_displacement_mpi   !displacement
          else
             gI_Indiv_Tasks_mpi%a(ipos(1),1) = gI_Indiv_Tasks_mpi%a(ipos(1)-1,2) + 1   !starting point
             if (ipos(2) <= gI_Vcolor%a(i1,4)) then
                gI_Indiv_Tasks_mpi%a(ipos(1),3) = gI_Vcolor%a(i1,3)       !count
             else
                gI_Indiv_Tasks_mpi%a(ipos(1),3) = gI_Vcolor%a(i1,3)-1       !count 
             endif
             gI_Indiv_Tasks_mpi%a(ipos(1),2) = gI_Indiv_Tasks_mpi%a(ipos(1),1) + gI_Indiv_Tasks_mpi%a(ipos(1),3) - 1  !end point
             gI_Indiv_Tasks_mpi%a(ipos(1),4) = gI_Indiv_Tasks_mpi%a(ipos(1)-1,4) + gI_Indiv_Tasks_mpi%a(ipos(1)-1,3)  !displacement 
          endif
       enddo !i2
    enddo !i1       
       
    if (new_id_mpi==6) then
       open(10,file=trim(cwd)//'output/test_IndivTasks.txt')
       do i1 = 1, gI_Indiv_Tasks_mpi%d1
          write(10,'(A,I10,A,4I8)') "gI_Indiv_Tasks_mpi%a(", i1, ',1:4) = ', gI_Indiv_Tasks_mpi%a(i1,:)
       enddo
       close(10)
    endif
    
    call sub_sgridA_set()  !set asset grid
       
    pinitA_domain(1,1:2) = (/sgridA%a(N_0-N_0+1,1), min(25.0d0, sgridA%a(N_0-N_0+1,VN_A(N_0)))/)
    
      
 END SUBROUTINE INIT_GRID

 SUBROUTINE INIT_GRID_IW() !initialize grids for IW decomposition
    IMPLICIT NONE
    integer(4) :: i1

       dvV_EV0%d1 = N_A
       dvV_EV0%d2 = N_AIME
       dvV_EV0%d3 = N_HEALTH
       dvV_EV0%d4 = N_RECSS
       dvV_EV0%d5 = N_FS
       dvV_EV0%d6 = N_H
       dvV_EV0%d7 = N_T_grid
       ALLOCATE( dvV_EV0%a( dvV_EV0%d1, dvV_EV0%d2, dvV_EV0%d3, dvV_EV0%d4, dvV_EV0%d5, dvV_EV0%d6, dvV_EV0%d7) )
       
       do i1 = 0, Ipt
           !policy functions    
           dcI0(i1)%d1 = N_A
           dcI0(i1)%d2 = N_AIME
           dcI0(i1)%d3 = N_HEALTH
           dcI0(i1)%d4 = N_RECSS
           dcI0(i1)%d5 = N_FS
           dcI0(i1)%d6 = N_H
           dcI0(i1)%d7 = dvV_EV0%d7
           ALLOCATE( dcI0(i1)%a( dcI0(i1)%d1, dcI0(i1)%d2, dcI0(i1)%d3, dcI0(i1)%d4, dcI0(i1)%d5, dcI0(i1)%d6, dcI0(i1)%d7) )
           
           !threshold value
           depsilon0(i1)%d1 = N_A
           depsilon0(i1)%d2 = N_AIME
           depsilon0(i1)%d3 = N_HEALTH
           depsilon0(i1)%d4 = N_RECSS
           depsilon0(i1)%d5 = N_FS
           depsilon0(i1)%d6 = N_H
           depsilon0(i1)%d7 = dvV_EV0%d7
           ALLOCATE( depsilon0(i1)%a( depsilon0(i1)%d1, depsilon0(i1)%d2, depsilon0(i1)%d3,  &
                   & depsilon0(i1)%d4, depsilon0(i1)%d5, depsilon0(i1)%d6, depsilon0(i1)%d7) )
       enddo !do i1
       
       do i1 = 0, Ipt+1
          dvV0(i1)%d1 = N_A
          dvV0(i1)%d2 = N_AIME
          dvV0(i1)%d3 = N_HEALTH
          dvV0(i1)%d4 = N_RECSS
          dvV0(i1)%d5 = N_FS
          dvV0(i1)%d6 = N_H
          dvV0(i1)%d7 = dvV_EV0%d7
          ALLOCATE( dvV0(i1)%a( dvV0(i1)%d1, dvV0(i1)%d2, dvV0(i1)%d3, dvV0(i1)%d4, dvV0(i1)%d5, dvV0(i1)%d6, dvV0(i1)%d7) )
       
          dcc0(i1)%d1 = N_A  !consumption policy function
          dcc0(i1)%d2 = N_AIME
          dcc0(i1)%d3 = N_HEALTH
          dcc0(i1)%d4 = N_RECSS
          dcc0(i1)%d5 = N_FS
          dcc0(i1)%d6 = N_H
          dcc0(i1)%d7 = dvV_EV0%d7
          ALLOCATE( dcc0(i1)%a( dcc0(i1)%d1, dcc0(i1)%d2, dcc0(i1)%d3, dcc0(i1)%d4, dcc0(i1)%d5, dcc0(i1)%d6, dcc0(i1)%d7) )
       
       enddo  !i1
       
       dvV_EV0%a(:,:,:,:,:,:,:) = dvV_EV%a(:,:,:,:,:,:,:)
       do i1 = 0, Ipt
           dcI0(i1)%a(:,:,:,:,:,:,:) = dcI(i1)%a(:,:,:,:,:,:,:)
           depsilon0(i1)%a(:,:,:,:,:,:,:) = depsilon(i1)%a(:,:,:,:,:,:,:)
       enddo
       do i1 = 0, Ipt+1
           dvV0(i1)%a(:,:,:,:,:,:,:) = dvV(i1)%a(:,:,:,:,:,:,:)
           dcc0(i1)%a(:,:,:,:,:,:,:) = dcc(i1)%a(:,:,:,:,:,:,:)
       enddo
    
  END SUBROUTINE INIT_GRID_IW     

  !The payroll taxes include the Social Security portion, 6.2% capped at $87, 900, 
  ! and the Medicare tax, 1.45% uncapped.
  SUBROUTINE INIT_TAX(id) !initialize TAX information
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: id
    real(8) :: rtax_ss, rtax_medicare, rtemp
    integer(4) :: i
    
    if ( (gcI_CF==3) .OR. (gcI_CF==7) .OR. (gcI_CF==10) .OR. (gcI_CF==73)) then
       rtax_ss = 0.0d0
       rtax_medicare = 0.0d0
    else
       rtax_ss = 0.062d0
       rtax_medicare = 0.0145d0 
    endif
    
    if (gcI_CF == 8) then
        rtemp = 1.5d0
    else
        rtemp = 1.0d0 
    endif
        
    do i = N_0, N_T        
        !gr_ftax_MTR: Federal tax, Marginal Tax Rate
        gr_ftax_MPTR(i,1) = rtax_ss+rtax_medicare
        gr_ftax_MPTR(i,2) = rtemp*0.1d0+rtax_ss+rtax_medicare
        gr_ftax_MPTR(i,3) = rtemp*0.15d0+rtax_ss+rtax_medicare
        gr_ftax_MPTR(i,4) = rtemp*0.25d0+rtax_ss+rtax_medicare
        gr_ftax_MPTR(i,5) = rtemp*0.25d0+rtax_medicare
        gr_ftax_MPTR(i,6) = rtemp*0.28d0+rtax_medicare
        gr_ftax_MPTR(i,7) = rtemp*0.33d0+rtax_medicare
        gr_ftax_MPTR(i,8) = rtemp*0.35d0+rtax_medicare
        gr_ftax_MPTR(i,:) = 1.0d0 - gr_ftax_MPTR(i,:)
    
        !gr_ftax_pbk: federal tax, bracket value
        gr_ftax_pbk(i,1) = 0.0d0  !Federal tax, including SS tax already
        gr_ftax_pbk(i,2) = 10250.0d0   !exemption+deduction = 3100.0d0+7150.0d0
        gr_ftax_pbk(i,3) = 20450.0d0
        gr_ftax_pbk(i,4) = 49150.0d0
        gr_ftax_pbk(i,5) = 87900.0d0    
        gr_ftax_pbk(i,6) = 110750.0d0
        gr_ftax_pbk(i,7) = 172950.0d0
        gr_ftax_pbk(i,8) = 329350.0d0
        gr_ftax_pbk(i,:) = gr_ftax_pbk(i,:)/grTotHours
    
        gr_ftax_pbase(i,1) = gr_ftax_pbk(i,1)
        gr_ftax_pbase(i,2) = gr_ftax_pbase(i,1)+(gr_ftax_pbk(i,2)-gr_ftax_pbk(i,1))*gr_ftax_MPTR(i,1)  !9465.88d0/grTotHours
        gr_ftax_pbase(i,3) = gr_ftax_pbase(i,2)+(gr_ftax_pbk(i,3)-gr_ftax_pbk(i,2))*gr_ftax_MPTR(i,2)  !17865.58d0/grTotHours 
        gr_ftax_pbase(i,4) = gr_ftax_pbase(i,3)+(gr_ftax_pbk(i,4)-gr_ftax_pbk(i,3))*gr_ftax_MPTR(i,3)  !40065.03d0/grTotHours 
        gr_ftax_pbase(i,5) = gr_ftax_pbase(i,4)+(gr_ftax_pbk(i,5)-gr_ftax_pbk(i,4))*gr_ftax_MPTR(i,4)  !66163.15d0/grTotHours 
        gr_ftax_pbase(i,6) = gr_ftax_pbase(i,5)+(gr_ftax_pbk(i,6)-gr_ftax_pbk(i,5))*gr_ftax_MPTR(i,5)  !82969.33d0/grTotHours 
        gr_ftax_pbase(i,7) = gr_ftax_pbase(i,6)+(gr_ftax_pbk(i,7)-gr_ftax_pbk(i,6))*gr_ftax_MPTR(i,6)  !126851.43d0/grTotHours 
        gr_ftax_pbase(i,8) = gr_ftax_pbase(i,7)+(gr_ftax_pbk(i,8)-gr_ftax_pbk(i,7))*gr_ftax_MPTR(i,7)  !229371.63d0/grTotHours 
        !This is to calculate the post-tax income directly
    enddo !do i
    
        
    !taxable social security benefit
    gr_sstax_1 = 2.5d4/grTotHours
    gr_sstax_2 = 9.0d3/grTotHours
    gr_sstax_3 = 3.4d4/grTotHours
    
    gr_sse_cap = 87900.0d0/grTotHours

  END SUBROUTINE INIT_TAX

  SUBROUTINE INIT_TAX_CF(t, rtax) !initialize TAX information, using modified tax rate rtax from t (inclusive) 
    IMPLICIT NONE
    INTEGER(4), INTENT(IN) :: t
    real(8), intent(in) :: rtax
    real(8) :: rtax_ss, rtax_medicare, rtemp
    integer(4) :: i
    
    rtax_ss = 0.062d0
    rtax_medicare = 0.0145d0 
    
    rtemp = 1.0d0 
    do i = N_0, t-1
        !gr_ftax_MTR: Federal tax, Marginal Tax Rate
        gr_ftax_MPTR(i,1) = rtax_ss+rtax_medicare
        gr_ftax_MPTR(i,2) = rtemp*0.1d0+rtax_ss+rtax_medicare
        gr_ftax_MPTR(i,3) = rtemp*0.15d0+rtax_ss+rtax_medicare
        gr_ftax_MPTR(i,4) = rtemp*0.25d0+rtax_ss+rtax_medicare
        gr_ftax_MPTR(i,5) = rtemp*0.25d0+rtax_medicare
        gr_ftax_MPTR(i,6) = rtemp*0.28d0+rtax_medicare
        gr_ftax_MPTR(i,7) = rtemp*0.33d0+rtax_medicare
        gr_ftax_MPTR(i,8) = rtemp*0.35d0+rtax_medicare
        gr_ftax_MPTR(i,:) = 1.0d0 - gr_ftax_MPTR(i,:)
    
        !gr_ftax_pbk: federal tax, bracket value
        gr_ftax_pbk(i,1) = 0.0d0  !Federal tax, including SS tax already
        gr_ftax_pbk(i,2) = 10250.0d0   !exemption+deduction = 3100.0d0+7150.0d0
        gr_ftax_pbk(i,3) = 20450.0d0
        gr_ftax_pbk(i,4) = 49150.0d0
        gr_ftax_pbk(i,5) = 87900.0d0    
        gr_ftax_pbk(i,6) = 110750.0d0
        gr_ftax_pbk(i,7) = 172950.0d0
        gr_ftax_pbk(i,8) = 329350.0d0
        gr_ftax_pbk(i,:) = gr_ftax_pbk(i,:)/grTotHours
    
        gr_ftax_pbase(i,1) = gr_ftax_pbk(i,1)
        gr_ftax_pbase(i,2) = gr_ftax_pbase(i,1)+(gr_ftax_pbk(i,2)-gr_ftax_pbk(i,1))*gr_ftax_MPTR(i,1)  !9465.88d0/grTotHours
        gr_ftax_pbase(i,3) = gr_ftax_pbase(i,2)+(gr_ftax_pbk(i,3)-gr_ftax_pbk(i,2))*gr_ftax_MPTR(i,2)  !17865.58d0/grTotHours 
        gr_ftax_pbase(i,4) = gr_ftax_pbase(i,3)+(gr_ftax_pbk(i,4)-gr_ftax_pbk(i,3))*gr_ftax_MPTR(i,3)  !40065.03d0/grTotHours 
        gr_ftax_pbase(i,5) = gr_ftax_pbase(i,4)+(gr_ftax_pbk(i,5)-gr_ftax_pbk(i,4))*gr_ftax_MPTR(i,4)  !66163.15d0/grTotHours 
        gr_ftax_pbase(i,6) = gr_ftax_pbase(i,5)+(gr_ftax_pbk(i,6)-gr_ftax_pbk(i,5))*gr_ftax_MPTR(i,5)  !82969.33d0/grTotHours 
        gr_ftax_pbase(i,7) = gr_ftax_pbase(i,6)+(gr_ftax_pbk(i,7)-gr_ftax_pbk(i,6))*gr_ftax_MPTR(i,6)  !126851.43d0/grTotHours 
        gr_ftax_pbase(i,8) = gr_ftax_pbase(i,7)+(gr_ftax_pbk(i,8)-gr_ftax_pbk(i,7))*gr_ftax_MPTR(i,7)  !229371.63d0/grTotHours 
        !This is to calculate the post-tax income directly
    enddo !do i
    
    rtemp = rtax    
    do i = t, N_T        
        !gr_ftax_MTR: Federal tax, Marginal Tax Rate
        gr_ftax_MPTR(i,1) = rtax_ss+rtax_medicare
        gr_ftax_MPTR(i,2) = rtemp*0.1d0+rtax_ss+rtax_medicare
        gr_ftax_MPTR(i,3) = rtemp*0.15d0+rtax_ss+rtax_medicare
        gr_ftax_MPTR(i,4) = rtemp*0.25d0+rtax_ss+rtax_medicare
        gr_ftax_MPTR(i,5) = rtemp*0.25d0+rtax_medicare
        gr_ftax_MPTR(i,6) = rtemp*0.28d0+rtax_medicare
        gr_ftax_MPTR(i,7) = rtemp*0.33d0+rtax_medicare
        gr_ftax_MPTR(i,8) = rtemp*0.35d0+rtax_medicare
        gr_ftax_MPTR(i,:) = 1.0d0 - gr_ftax_MPTR(i,:)
    
        !gr_ftax_pbk: federal tax, bracket value
        gr_ftax_pbk(i,1) = 0.0d0  !Federal tax, including SS tax already
        gr_ftax_pbk(i,2) = 10250.0d0   !exemption+deduction = 3100.0d0+7150.0d0
        gr_ftax_pbk(i,3) = 20450.0d0
        gr_ftax_pbk(i,4) = 49150.0d0
        gr_ftax_pbk(i,5) = 87900.0d0    
        gr_ftax_pbk(i,6) = 110750.0d0
        gr_ftax_pbk(i,7) = 172950.0d0
        gr_ftax_pbk(i,8) = 329350.0d0
        gr_ftax_pbk(i,:) = gr_ftax_pbk(i,:)/grTotHours
    
        gr_ftax_pbase(i,1) = gr_ftax_pbk(i,1)
        gr_ftax_pbase(i,2) = gr_ftax_pbase(i,1)+(gr_ftax_pbk(i,2)-gr_ftax_pbk(i,1))*gr_ftax_MPTR(i,1)  !9465.88d0/grTotHours
        gr_ftax_pbase(i,3) = gr_ftax_pbase(i,2)+(gr_ftax_pbk(i,3)-gr_ftax_pbk(i,2))*gr_ftax_MPTR(i,2)  !17865.58d0/grTotHours 
        gr_ftax_pbase(i,4) = gr_ftax_pbase(i,3)+(gr_ftax_pbk(i,4)-gr_ftax_pbk(i,3))*gr_ftax_MPTR(i,3)  !40065.03d0/grTotHours 
        gr_ftax_pbase(i,5) = gr_ftax_pbase(i,4)+(gr_ftax_pbk(i,5)-gr_ftax_pbk(i,4))*gr_ftax_MPTR(i,4)  !66163.15d0/grTotHours 
        gr_ftax_pbase(i,6) = gr_ftax_pbase(i,5)+(gr_ftax_pbk(i,6)-gr_ftax_pbk(i,5))*gr_ftax_MPTR(i,5)  !82969.33d0/grTotHours 
        gr_ftax_pbase(i,7) = gr_ftax_pbase(i,6)+(gr_ftax_pbk(i,7)-gr_ftax_pbk(i,6))*gr_ftax_MPTR(i,6)  !126851.43d0/grTotHours
        gr_ftax_pbase(i,8) = gr_ftax_pbase(i,7)+(gr_ftax_pbk(i,8)-gr_ftax_pbk(i,7))*gr_ftax_MPTR(i,7)  !229371.63d0/grTotHours
        !This is to calculate the post-tax income directly
    enddo !do i


  END SUBROUTINE INIT_TAX_CF

  
 !initialize HEALTH information. call after INIT_GRID
 !called in the master, head slaves and slaves
 SUBROUTINE INIT_HEALTH(id) 
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: id
    integer(4) :: i
    
    if (IHasHealthAge == 0) then !in the baseline model without health status
        MSM_Indiv_Ihealth%a(N_0-N_0+1:,1:N_Sim) = 1  !there is only ONE health status
        return
    endif
    
       do i = 1, N_Sim
                    
          !health status is Good (2) from N_0 to IHasHealthAge-1.  
          MSM_Indiv_Ihealth%a(N_0-N_0+1:IHasHealthAge-1-N_0+1,i) = 2
          ! At IHasHealthAge, the health is randomly determined.
          
          ! to determine health status at IHasHealthAge
          git = IHasHealthAge - 1
          git_pos = git - N_0 + 1
          if (MSM_Indiv%a(11,git_pos,i) <= gr_health_matrix(git,1,0)) then
             MSM_Indiv_Ihealth%a(git_pos+1,i) = 1  !health is excellent
          elseif ( MSM_Indiv%a(11,git_pos,i) <= sum(gr_health_matrix(git,1:2,0)) ) then
             MSM_Indiv_Ihealth%a(git_pos+1,i) = 2  !health is good
          elseif ( MSM_Indiv%a(11,git_pos,i) <= sum(gr_health_matrix(git,1:3,0)) ) then
             MSM_Indiv_Ihealth%a(git_pos+1,i) = 3  !health is bad
          else   
             MSM_Indiv_Ihealth%a(git_pos+1,i) = 4  !health is disable
          endif
                
          ! to determine health status from IHasHealthAge+1 to N_T
          do git = IHasHealthAge, N_T-1
             git_pos = git - N_0 + 1 
             giHEALTH = MSM_Indiv_Ihealth%a(git_pos,i)  !health status, 1-Excellent, 2-Good, 3-bad, 4-Disab 
             if (giHEALTH<4) then                  !being in non-disable (excellent 1, good 2, or bad 3) health
                if (MSM_Indiv%a(11,git_pos,i) <= gr_health_matrix(git,giHEALTH,1)) then
                   MSM_Indiv_Ihealth%a(git_pos+1,i) = 1  !health becomes excellent
                elseif ( MSM_Indiv%a(11,git_pos,i) <= sum(gr_health_matrix(git,giHEALTH,1:2)) ) then
                   MSM_Indiv_Ihealth%a(git_pos+1,i) = 2  !health becomes good
                elseif ( MSM_Indiv%a(11,git_pos,i) <= sum(gr_health_matrix(git,giHEALTH,1:3)) ) then
                   MSM_Indiv_Ihealth%a(git_pos+1,i) = 3  !health becomes bad
                else   
                   MSM_Indiv_Ihealth%a(git_pos+1,i) = 4  !health becomes disable
                endif    
             else  !being in disable (4) health
                MSM_Indiv_Ihealth%a(git_pos+1,i) = 4   !disable is an absorbing status
             endif !if (giHEALTH<4) then
          enddo !do git
          
          !now set health status according to gcI_CF
          if ( (gcI_CF > 40) .AND. (gcI_CF < 60)) then !41-45, 46-50, 51-55: anticipated or un-anticipated
             select case (gcI_CF)
                case (41,46,51)
                   MSM_Indiv_Ihealth%a(51-N_0+1:,i) = MSM_Indiv_Ihealth%a(50-N_0+1,i)  !health stays unchanged after 50
                case (42,47,52)
                   MSM_Indiv_Ihealth%a(50-N_0+1:,i) = 1 !health becomes excellent and stays since 50
                case (43,48,53)
                   MSM_Indiv_Ihealth%a(50-N_0+1:,i) = 2 !health becomes good and stays good since 50
                case (44,49,54)
                   MSM_Indiv_Ihealth%a(50-N_0+1:,i) = 3 !health becomes bad and stays bad since 50
                case (45,50,55)
                   MSM_Indiv_Ihealth%a(50-N_0+1:,i) = 4 !health becomes disabled and stays disabled since 50   
             end select 
          endif  !if (gcI_CF > 0)                        
       enddo !do i = 1, N_Sim
    
 END SUBROUTINE INIT_HEALTH
 
 !initialize Family Status information. call after INIT_GRID
 !called in the master, head slaves and slaves
 !family status, 1-single, 2-married, spouse not working; 3-married, spouse working
 SUBROUTINE INIT_FamilyStatus(id) 
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: id
    integer(4) :: i, j
    real(8) :: rtemp
    
    do i = 1, N_Sim
       if (gcL_noFS) then
          MSM_Indiv_IFS%a(:,i) = 1  !single always
          MSM_Indiv%a(16,:,i) = 0.0d0  !spousal income
          cycle
       endif !if (gcL_noFS)
           
       MSM_Indiv_IFS%a(N_0-N_0+1,i) = 1  !single at N_0
       !Family status next period
       do git = N_0, N_T-1
          git_pos = git - N_0 + 1 
          giFS = MSM_Indiv_IFS%a(git_pos,i)
          
          rtemp = 0.0d0
          do j = 1, N_FS
             rtemp = rtemp + gr_FS_trans(git,giFS,j)
             if (MSM_Indiv%a(15,git_pos,i) <= rtemp) then
                MSM_Indiv_IFS%a(git_pos+1,i) = j
                exit
             endif
          enddo !do j
          if (j > N_FS) MSM_Indiv_IFS%a(git_pos+1,i) = N_FS
       enddo !do git

       !spouse income
       do git = N_0, N_T
          git_pos = git - N_0 + 1 
          if (MSM_Indiv_IFS%a(git_pos,i) == 3) then   !spouse working
             if (gcL_spinc_discrete) then
                rtemp = 0.0d0 
                do j = 1, N_FS_spinc_int
                   rtemp = rtemp + gr_FS_spinc_xw(git,j,2)
                   if (MSM_Indiv%a(16,git_pos,i) <= rtemp) then
                      MSM_Indiv%a(16,git_pos,i) = gr_FS_spinc_xw(git,j,3)
                      exit
                   endif
                   if (j>N_FS_spinc_int) MSM_Indiv%a(16,git_pos,i) = gr_FS_spinc_xw(git,N_FS_spinc_int,3)
                enddo !do i2 
             else
                rtemp = gr_FS_inc(git,1) + gr_FS_inc(git,2)*MSM_Indiv%a(16,git_pos,i) 
                MSM_Indiv%a(16,git_pos,i) = min(exp(rtemp), 2.0d0*maxval(gr_FS_spinc_xw(git,:,3)))         !ln(spinc)
             endif  !if (gcL_spinc_discrete)
          else
             MSM_Indiv%a(16,git_pos,i) = 0.0d0
          endif
       enddo  !do git
    enddo !do i = 1, N_Sim
 END SUBROUTINE INIT_FamilyStatus
  
  !value/policy functions and some values initilization, in all nodes
  !This is called in each iteration
  subroutine sub_node_init()
      implicit none
      integer(4) :: i 
      real(8) :: rtemp
      
      !!!! !set environment
      if (gcI_CF .NE. 0) then
         gcL_IsOutput = .TRUE.
         gcL_cal_shadowprice = .FALSE.
         pw = 1.0d0
         gcI_elas_t = 0
         gcI_taxhike_t = 0
         pNRA = 65  !restore NRA to 65
         
         pPIAratio(1:3) = (/0.9d0, 0.32d0, 0.15d0/)  ! PIA raito in Equation C4: 0.9, 0.32, 0.15 in the data
    
         select case (gcI_CF)
            case (2, 9)               
               pNRA = 67
            case (-10)   !anticipated elasticity calculation
               gcI_elas_t = NINT(gr_Vpar_mpi(N_PAR+2))
               pw(gcI_elas_t) = 1.1d0 
            case (-12)   !unanticipated elasticity calculation
               gcI_elas_t = NINT(gr_Vpar_mpi(N_PAR+2))
            case (-15)   !tax hike, anticipated
               gcI_taxhike_t = NINT(gr_Vpar_mpi(N_PAR+2)) 
               rtemp = 1.5d0
               call INIT_TAX_CF(gcI_taxhike_t, rtemp)
            case (20)
               pPIAratio(1:3) = (/0.61d0, 0.61d0, 0.15d0/)
            case (21)
               pPIAratio(1:3) = 0.1d0
            case (22)
               pPIAratio(1:3) = 0.3d0
            case (23)
               pPIAratio(1:3) = 0.5d0
            case (24)
               pPIAratio(1:3) = 0.7d0
            case (25)
               pPIAratio(1:3) = 0.9d0
            case (28)
                pPIAratio(1:3) = 0.4475d0
            case (70, 88)
            case default
         end select
         if (gcI_CF .NE. -15) then
             call INIT_TAX(myid_mpi)
         endif
         
 
         !initialize health status for simulation
         !gcI_CF: 0 not, 41-45, 46-50, 51-55
         !ATTENTION: these experiments should be executed in the very end.
         if (IHasHealthAge>=1 .AND. gcI_CF>=41 .AND. gcI_CF<=55) call INIT_HEALTH(myid_mpi) 
      endif !if (gcI_CF .NE. 0)
      
      gcI_itercount = gcI_itercount + 1  !iteration counting in head slaves and slaves
       
      !!!! !initialize policy and value
      
          do i = 0, Ipt+1 
              dcc(i)%a = 0.0d0
              dvV(i)%a = - NR_BIG2
              dciRECSS(i)%a = 0
          enddo
          dvV_EV%a = - NR_BIG2
          do i = 0, Ipt
              dcI(i)%a = 0.0d0
              depsilon(i)%a = 0.0d0
          enddo
  end subroutine sub_node_init
 
 

SUBROUTINE DESTROY_GRID(id)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: id
    integer(4) :: i 
    
    if (id == master_mpi) then
       DEALLOCATE(MSM_Moments_Data%a)
       DEALLOCATE(MSM_wgt%a)
       DEALLOCATE(gmmse_S%a)
       deallocate(MSM_wgt_scale%a)
    endif
    
    !all nodes, including master and slaves
    deallocate(gI_VTasks_mpi%a) 


    
    if (gcL_debugprint) write(*,'(A,I4,A,I4,A,I4,A)') "iter#", gcI_itercount, ' slave #', new_id_mpi, &
                           & ' in color# ',color_mpi, ' destroying 1'
    DEALLOCATE(dvV_EV%a)    
    if (gcL_debugprint) write(*,'(A,I4,A,I4,A,I4,A)') "iter#", gcI_itercount, ' slave #', new_id_mpi, &
                           & ' in color# ',color_mpi, ' destroying 3'
   
    do i = 0, Ipt
       DEALLOCATE(dcI(i)%a)
       DEALLOCATE(depsilon(i)%a)
    enddo

    do i = 0, Ipt+1
       DEALLOCATE(dvV(i)%a)
       DEALLOCATE(dcc(i)%a)
       DEALLOCATE(dciRECSS(i)%a)
    enddo
    
    deallocate(MSM_sim4m_i%a)
    DEALLOCATE(MSM_Moments_Sim%a)
    
       
    DEALLOCATE(MSM_Indiv_initX%a)
    DEALLOCATE(MSM_Indiv%a)    
    DEALLOCATE(MSM_Indiv_Ihet%a)
    DEALLOCATE(MSM_Indiv_Irecss%a)   
    DEALLOCATE(MSM_Indiv_Ihealth%a)
    DEALLOCATE(MSM_Indiv_IFS%a)
    
    deallocate(gI_Vcolor%a)
    
    DEALLOCATE(sgridA%a)
    DEALLOCATE(sgridAnrlb%a)
    DEALLOCATE(gr_Amax%a)
    DEALLOCATE(sgridH%a)
    DEALLOCATE(sgridAIME%a)
    DEALLOCATE(sgridRECSS%a)
    
    if (gcL_searchBC .OR. (N_pbequest == 2)) then
        DEALLOCATE(sgridA_BC%a)
        DEALLOCATE(sgridAnrlb_BC%a)
    endif
    
END SUBROUTINE DESTROY_GRID

SUBROUTINE DESTROY_GRID_IW()
    IMPLICIT NONE
    integer(4) :: i
    
    DEALLOCATE(dvV_EV0%a)
    
    do i = 0, Ipt
       DEALLOCATE(dcI0(i)%a)
       DEALLOCATE(depsilon0(i)%a)
    enddo

    do i = 0, Ipt+1
       DEALLOCATE(dvV0(i)%a)
       DEALLOCATE(dcc0(i)%a)
    enddo    
END SUBROUTINE DESTROY_GRID_IW

SUBROUTINE DESTROY_VALUE_DATA(id)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: id
    
    deallocate(cps_aimeshare%a) 
        
END SUBROUTINE DESTROY_VALUE_DATA

END MODULE INIT


